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Abstract 

A one-way street of width M is modeled as a set of M parallel one- 
dimensional TASEPs. The intersection of two perpendicular streets 
is a square lattice of size M x M. We consider hard core particles 
entering each street with an injection probability a. On the inter- 
section square the hard core exclusion creates a many-body problem 
of strongly interacting TASEPs and we study the collective dynamics 
that arises. We construct an efficient algorithm that allows for the sim- 
ulation of streets of infinite length, which have sharply defined critical 
jamming points. The algorithm employs the 'frozen shuffle update', 
in which the randomly arriving particles have fully deterministic bulk 
dynamics. High precision simulations for street widths up to M = 24 
show that when a increases, there occur jamming transitions at a se- 
quence of M critical values a 1 ^ < a 1 ^ I I _ 1 < . . . < a^ 1 . As M grows, 
the principal transition point decreases roughly as ~ (logM) -1 in 
the range of M values studied. We show that a suitable order param- 
eter is provided by a reflection coefficient associated with the particle 
current in each TASEP. 

Keywords: exclusion process, crossing flows, pedestrian traffic, frozen 
shuffle update 
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1 Introduction 



Pedestrian motion in dense environments is of both theoretical and 
practical interest. Instances of applications are shopping streets, 
waiting lines, crowds that enter or leave a confined space, and so 
on. Under such circumstances simplified models may help understand the 
behavior of individuals as well as the collective behavior that results from it 
[H 121 [31 H] . A particular class of such models is based on cellular automata 
[HI El IE] • For modeling unidirectional one- dimensional traffic, whether 
it be of particles, vehicles, or pedestrians [9], one popular tool is the To- 
tally Asymmetric Simple Exclusion Process (TASEP) [TO]. This stochastic 
process belongs to the larger class of models of random walkers with hard 
core interactions and various different TASEP versions have been studied 
by physicists and mathematicians alike since several decades. Specifically, a 
TASEP is a system of hard core particles that advance along a linear lat- 
tice in a single direction. The TASEP may be used as a building block for 
traffic flow models in more complicated geometries. There is a rich litera- 
ture on junctions and bifurcations [HI [121 EES EH]; as we U as on intersections 

[IS [m El EH [HI [20] of TASEPs - In those a PP roaches the interaction be- 
tween the different TASEPs is 'weak' in the sense that they are coupled only 

on a set of generally well-separated sites. Strongly interacting TASEPs were 

considered, in particular, in the so-called BHL model proposed by Biham et 

al. [21], which in turn has given rise to an offspring of variants on some of 

which we will comment. 

Although pedestrian traffic flow is our basic motivation, we present below 
a model which, because of its simplicity, has an intrinsic interest that extends 
beyond this particular application. From a wider perspective it is an example 
of a driven nonequilibrium system, and, contrary to many other such systems 
that have been studied in the literature and that have stochastic dynamics 
[2"2"l |2"3] , this one is deterministic. 

We define a street of width M as a set of M parallel linear lattices to 
be called lanes, each of which carries a TASEP. We wish to study what 
happens when two such streets, both thought of as being of infinite length, 
intersect perpendicularly. Figure [1] shows the geometry for the finite length 
case. Particles in the horizontal street are injected from the left onto empty 
sites with a probability a per time step. Particles in the vertical street 
are injected from below with the same probability. The resulting average 
incoming current in a lane will be denoted by J free (a). 

The intersection area of the two streets is a finite M x M square lat- 
tice. On this square we have a problem of strongly interacting TASEPs, 
whose study is the purpose of this work. Lanes are numbered by an index 
m = 1, 2, . . . , M from the outer ones inward. Unless symmetry were bro- 
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ken spontaneously, the mth horizontal and mth vertical lane are statistically 
identical. One basic question here, as in any other traffic flow model, is to 
determine the outgoing currents J^{oi) as a function of a when the system 
is in a stationary state. Contrary to the incoming current, the outflow must 
be expected to be m and M dependent. A lane is said to be in a free flow 
phase when its outgoing current is equal to the incoming current, that is, 
when J^f(a) = J frcc (a). At sufficiently low a free flow is to be anticipated 
for all m; however, it is inevitable that when a increases above a certain 
threshold, part or all of the lanes undergo a jamming transition. 

The update scheme is an essential part of the definition of any TASEP. 
We choose to employ here the frozen shuffle update [2U [25], [20], whose char- 
acteristic is that in each time step the particles are updated according to a 
sequence fixed once and for all. Particles entering the system are inserted in 
this sequence and particles leaving are deleted from it. This update, orig- 
inally proposed as an alternative to parallel updat^E has the advantage of 
eliminating the algorithmic conflicts that arise when in the same time step 
two particles have the same target site. It conserves the advantage of the 
parallel update in that it bounds the fluctuations that random sequential 
update would cause. In our present implementation of the frozen update 
scheme all allowed moves are carried out with probability one and hence the 
bulk dynamics is that of a cellular automaton: it is fully deterministic. 

Under frozen shuffle update, once a particle has left the intersection square 
it continues unimpeded at unit speed, and hence the street segments beyond 
the intersection square need not be considered. The street segments leading 
up to this square, however, play an essential role, since they are the place 
where waiting lines may develop, whether temporarily or permanently. Such 
waiting lines are an intrinsic part of the traffic flow problem and we wish to 
fully account for them. In any simulation the ingoing street segments are 
necessarily of finite length L, as shown in figure [TJ A key observation is, 
however, that sharply defined jamming transition points on the a axis can 
exist only in the limit L — > oo. Indeed, as previous work [2U] has shown, 
there is a sharp transition point for crossing streets even of width M = 1, as 
long as their length L is infinite. This being so, the present work does the 
following. 

First, by combining theoretical and algorithmic arguments, we show that 
it is possible to integrate out the degrees of freedom in the half-way infinite 
street segments leading up to the intersection square. There then results an 
efficient algorithm that simulates the particle motion on the finite M x M 
square, subject to boundaries on its left and lower edge that represent the 

frozen shuffle update may also be seen as a variant of random shuffle update [261 127) ; 
in the latter a new random particle order is drawn before each time step. 
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incoming street segments of infinite length, L = oo. Finite size effects have 
thus been eliminated. The appropriate boundary conditions are derived and 
formulated in terms of memory variables to be defined in section [31 

Secondly, we present the results of the simulation of this interacting street 
model as a function of its two parameters a and M. We obtain the outgoing 
currents J^(a) and find that when a increases from to 1, there appear M 
critical values 

a%<*%_ x <...<< (1.1) 

at which, successively, the lanes with indices m = M, M — 1, . . . , 1 become 
jammed. We will refer to as the principal critical point of the size M 
intersection. For a beyond the critical value a^f, the outflow of particles 
in the mth horizontal and vertical lane cannot keep up with the incoming 
current and in those lanes ever growing waiting lines develop. 
We relate the incoming and outgoing current in lane m by 

J% = (1 " f^V free - (1-2) 

The product R^f J free in (11.21) may be interpreted as the reflected currentH in 
the mth lane and we therefore call R^f its reflection coefficient. We will show 
that there exists a simple theoretical relation between this coefficient and the 
memory variables occurring in the boundary condition. In our Monte Carlo 
work we deduce the critical point values from the a dependence of R^f 
through the criterion 

R£(«) = o, 

R^(a)>0, a>«„, m = l,2,...,M. (1.3) 

For a% < a < 1 the coefficient R^f (a) increases monotonously and we may 
appropriately consider it as an order parameter. The main results of our 
simulations are the determination of the critical points a% and of the curves 

R£(«). 

This paper is set up as follows. In section [2] we define the rules of motion 
for what we will call the 'full' algorithm, that is, the one executed on the 
lattice of figure [Tj, which has a finite L. In section [3] the memory variables 
are introduced; we then show how for L = oo the full algorithm gives rise to 
a 'reduced' algorithm, restricted to the M x M interaction square, and with 
boundary conditions formulated in terms of these new variables. In section 
H] we establish the theoretical expressions, valid for each lane separately, 
that relate the memory variable to the transmitted current and the speed 
of propagation of the reflected current. In section Owe present and discuss 
simulation results for lattices up to M = 24. In section [6] we make several 
final remarks and conclude. 

2 This even though all particles move only in a single direction; see section POl 
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Figure 1: Intersection of two one-way streets of width M. The x particles, moving 
to right, and the y particles, moving upward, are represented by blue and red 
triangles, respectively. The control parameter a is the injection probability per 
time step for each lane. Small arrows point towards the 'injection' sites. Lanes 
are numbered m = 1, 2, . . . , M. The region bordered by the heavy solid line is the 
'intersection square'. The 'entrance sites' defined in section [3] are shaded. 
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2 Intersecting street model 



We consider the lattice of figure [I] showing a 'horizontal' and a 'vertical' one- 
way street labeled 'x' and 'y', respectively, and each composed of M parallel 
lanes. A lane is identified by an index pair (m, 8) with m = 1,2, ...,M 
and 8 = 'x' or l y\ The length L of the street segments leading up to the 
intersection is a large finite numbeiE). The two streets intersect according 
to a square of M x M sites. Two types of hard-core particles move on the 
lattice, x particles arriving from the left along the horizontal street, and y 
particles arriving from below along the vertical street. New particles may 
be injected onto the sites of the leftmost column of the horizontal street 
and of the bottom row of the vertical street; we will refer to these as the 
injection sites. We let a > denote the injection rate of particles onto an 
empty injection site. The ensuing injection probability a = 1 — e~ a is the 
probability that an injection site, when empty at a certain instant of time t, 
will get occupied during the unit time interval (t, t + 1) that follows. We will 
employ a rather than a as the control parameter. 

2.1 Full algorithm 

The rules for the particle motion on this lattice are as follows. The algo- 
rithm is a succession of time steps labeled by an integer time step index 
s = 1, 2, 3, . . .. The sth time step takes the system from time t — s — 1 to 
t = s. It executes all events (particle moves, injections and deletions) that 
occur during the time interval s — 1 < t < s. 

With the frozen shuffle update [2U 123 ED] employed here, either at the 
initial time t = 0, or else when it enters the system, the ith particle in 
lane (m, 8) is assigned, in a way to be described below, a random phase 
T mS £ (0, 1); it keeps this phase as a fixed attribute. The phases define 
in which order the particles are updated at each time step. For practical 
reasons the algorithm decouples each time step update into a 'bulk' step and 
an 'injection' step, both to be described now. 

2.1.1 Bulk step 

The sth bulk time step consists of a sweep through the 'update sequence', 
that is a list, common to all lanes, of all particles present in the system at the 
beginning of that time step, ordered according to increasing phases. Thus, 
during the bulk step, each particle present in the system is visited exactly 
once. This bulk step parallels the evolution of the continuous time variable 
t, the turn for particle % to be updated being identified with the instant of 
time t = s - 1 + t™ 5 . 

3 In section |3] we will set L = oo. 
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General case. The case of a general particle % in a lane (m, S) is simple. 
When at time t = s — 1 + t™ s its turn to be updated has come, then if at 
that instant of time its target site is occupied, it does not move; if its target 
site is empty, it moves. 

First special case. A particle % occupying the last site of its lane (an 'exit 
site') has no target site. When its turn to be updated has come, it is deleted 
from the system; in particular, its phase is deleted from the update sequence. 
We will comment further on this boundary condition in section 12.21 

Second special case. When a particle j moves off the injection site of a 
lane (m, 5), say at time s — 1 + rj™ 5 , the algorithm needs to do more work. In 
that case, in order to prepare for the injection of the next particle, j + 1, onto 
that injection site, a random time interval TJ^ is drawn from the exponential 
distribution 

P(T) = ae- aT , T>0, (2.1) 
and j + 1 is scheduled to be injected at time tj+i given bjfl 

tj+1 = S - 1+Tj + T ™{ . (2.2) 

This implies, first, that this new particle has a phase r^+i equal to the frac- 
tional part of its time of injection, that is, 

r^^ + T™*) modi; (2.3) 

and, secondly, that its injection will occur during time step s J+ i = [tj+ij + 
1 = s + \jj + X^iJ , where [^J denotes the largest integer contained in x. 
Particle j + 1, its phase ry+i, its time step of injection Sj + i, and its lane index 
(m, S), are placed on a waiting list to which the algorithm will return during 
the injection step. This completes the discussion of the second special case. 

Remark (i). We have Sj+\ > s, that is, the injection may occur either 
during the same time step s in which the injection site was emptied, or during 
any later time step. 

Remark (ii). The injection rule implements independent and uniformly 
distributed arrival times subject to interparticle hard core exclusion [25| 120]. 

2.1.2 Injection step 

The waiting lists of each of the 2M lanes (m, 5) are consulted independently. 
If the particle, say k, on that list has a time step of injection Sk equal to the 
current time step s, then the algorithm implements the injection; the phase 
t™ 6 is inserted in the update sequence; and the injection is considered to 

4 When no confusion can arise we will suppress the superscript 'mJ' on t" iS . Similarly, 
we will not append this superscript on quantities like tj+i and fcj+i introduced in this 
subsection. 
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have occurred physically at time s — 1 + t™ 5 . In all following time steps, and 
until particle k has left the system, its position will be updated as part of 
the bulk step. 

Remark. In the absence of any obstacles to the particle flow the above 
injection procedure brings the incoming lane segment into a 'free flow con- 
figuration', that is, one in which each particle advances by one lattice unit 
each time stejfi This configuration has an average particle density p free (a) 
and an average current J frec (a) = vp (a) given by [251 120] 

/ rc » = T^, J frc » = T ^, (2.4) 
1 + a 1 + a 

where we used that v — 1. 



2.1.3 Initial state 

At time t = we initialize the system in a random free flow configuration in 
which the flow has just reached the entrance sites (shown shaded in figured]) 
of the intersection square, but has not entered the square itself yet. Such a 
configuration is conveniently generated as follows. 

(1) Set time equal to t — —L + 1. 

(2) On each of the injection sites of the 2M lanes independently, deposit 
a particle with the probability p iree (a) given by (12 .4p and leave it unoccupied 
with the complementary probability 1 — p (a). If a particle is deposited, 
assign to it a phase r™ 5 drawn randomly and uniformly from (0, 1). If no par- 
ticle is deposited, draw a T™ 5 from (12. ip and prepare to occupy the injection 
site at a time t = —L + 1 + T 1 m<5 (that is, during some time step s > —L + 2) 
by a particle of phase r] n<5 = T™ 5 mod 1. 

(3) Execute L — 1 time steps to take the system from time t = —L + 1 to 
t = 0. During these time steps the 2M lanes do not interact and no blocking 
occurs. As a result, each particle initially deposited on an injection site will at 
t = occupy the corresponding entrance site, and the lane segments leading 
up to the intersection square will carry a free flow configuration. 

Remark. Both the random initial state and the time evolution of the 
system are fully determined by the sequences of time intervals 

Tr 5 ,T 2 m<5 ,T 3 m5 ,..., m = l,2,...,M, 5 = V or V , (2.5) 

which, in particular, determine the t™ 5 for j = 1, 2, 3 . . . through (I2.3p . 

5 A 'free flow configuration' is a microstatc. Note that the macroscopic 'free flow phase' 
was defined in section [T] by the equality of inflow and outflow. 
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Figure 2: The M x M intersection square on which the 'reduced' algorithm is 
defined. An extra row and column of 'entrance sites' (shaded) along the left and 
lower border serve to take care of the boundary conditions. 

2.2 The exit boundary 

The lane segments beyond the intersection square are initially empty. The 
first particle to enter onto the first site of such a lane segment will from 
that moment on proceed unblocked at speed v = 1. Each next particle will 
enter the same segment only after its predecessor has left the first site and 
from then on similarly proceeds at speed v = 1. Hence each outgoing lane 
segment carries a free flow configuration and does not exert any feedback on 
the intersection square. The boundary condition that we applied in section 
12.1.11 - namely to eliminate a particle which when its update time has come 
no longer has a target site - therefore correctly represents an outgoing free 
flow heading to infinity. We will henceforth apply this boundary condition 
at the immediate upper and right border of the intersection square. 

3 The intersection square 

We construct in this section an algorithm for simulating the crossing street 
model which is an exact consequence of the above 'full' algorithm but which 
is much more efficient. This 'reduced' algorithm retains only the particle 
motion on the intersection square. As we have seen, the particles exiting this 
square cease to play a role. We will show below that it is possible - in the 
limit L — > oo (see figure [1]), which is precisely the one of interest - also to 
eliminate the incoming flows from the simulation. This is achieved with the 
aid of appropriate boundary conditions defined on the entrance sites shown 
shaded in figure [U Hence the system geometry becomes that of figure [2] 
For L = oo all finite size effects are eliminated and the system has sharply 
denned critical points. Any remaining uncertainty on the location of these 
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points will be due only to the finite duration of the simulation. 
3.1 Memory variables I mS (s) 

All statistical properties of the free flow injected into a lane (m, 5) are known. 
If this flow did not encounter any obstacles, then it would arrive unmodified 
at the entrance site L time steps later, and the boundary condition at that 
site would be known in a simple way. However, blocking of particles in the 
intersection square will typically create at the entrance site a waiting line 
of variable length which complicates the boundary conditions. This waiting 
line may be either ever growing or, instead, only fluctuating with a finite 
localization length, two possibilities that correspond to the lane being in the 
'jammed flow' and the 'free flow' phase, respectively. The discussion of the 
present section applies to both cases. 

The solution to the boundary condition problem will be shown to consist 
in introducing an auxiliary time dependent memory variable I m5 (s). We will 
be able by means of this variable to determine the arrival time t = tj of each 
particle j at the entrance site exactly as if, starting at time t — — oo, it had 
made its way through a street segment of length L = oo, taking correctly 
into account any time it may have spent blocked in the waiting line. 

We associate with each entrance site a 'reference particle': this will be 
the particle located on that site if there is one, and the first particle going 
to arrive on that site otherwise. Consider now the incoming segment of lane 
(m,5). During each time step every particle either advances by one lattice 
unit or incurs a 'unit time delay'. We define the memory variable I mS (s) 
as the accumulated time delay that the reference particle in lane (m, S) has 
incurred in the course of its history up to and including the sth time step. 
It should be noted that I mS (s) is not associated with a fixed particle; when 
a particle leaves the entrance site, the next one in the same lane takes over 
its role as the reference particle and continues to carry the variable I mS (s), 
the value of which is in general affected by this takeover process. The I m5 (s) 
are nonnegative integers; for the initial state described in subsection 12. 1 . 31 we 
have I mS (0) = 0. 

In order to determine the time evolution of the 7 m<5 (s) we examine in 
detail how the sth time step relates I mS (s) to I mS (s — 1). In any time step 
three different events may occur at the entrance site of lane (m, 5); we denote 
them by the symbols 'B' (for 'blocked'), 'A' (for 'advancing'), and 'E' (for 
'empty'). 

Event B. The reference particle occupies the entrance site but is blocked 
and does not move. In that case its delay I mS (s) is augmented by one unit, 

I mS ( s ) = I mS ( s - 1) + 1. (3.1a) 
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Event A. The reference particle, say j, occupies the entrance site and 
during the sth time step advances into the intersection square. The next 
particle, j + 1, becomes the reference particle and J m<5 (s), now redefined as 
the accumulated time delay incurred by j + 1, must be recalculated. If j 
would never have been blocked, its successor j + 1 would have performed the 
same jumps as j but with a time delay Tj^ + 1. Here TJ^ is the difference 
between the departure time of j and the arrival time of j + 1 on an arbitrary 
site, and the extra +1 is the sojourn time of an unblocked particle on that 
site. 

The extra headway T™^ 1 allows j to be blocked a total number of L^j+iJ 
times without blocking j + 1. Hence if the time delay I m5 (s — 1) of j is small 
enough to satisfy I mS (s — 1) < L^J+iJ) then I mS (s) = 0. Each supplementary 
blocking of j leads to a blocking, and hence to a unit time delay, of j + 1, 
and so if I m5 (s - 1) > [T™{ + lj, then I mS (s) = I m5 (s - 1) - [T™\\. This 
may be combined into 

I m5 (s) = max (I m5 (s - 1) - [T™{\ , 0). (3.1b) 

We note that in event A particle j + 1 may or may not arrive on the entrance 
site during the time step s under consideration. 

Event E. The reference particle does not occupy the entrance site at the 
beginning of the sth time step. Being free to move, it comes one lattice site 
closer to the entrance site during that time step. Therefore its incurred time 
delay remains unchanged, that is, 

I mS (s) = I m5 (s - 1). (3.1c) 

Again, in this event the reference particle may or may not arrive on the 
entrance site during time step s. 

Equations fl3.1ap -( )3.1cp govern the time evolution of I mS (s). We are led 
to the important conclusion that they involve exclusively the local motion of 
a single reference particle on or near the entrance site of lane (m, 5); there is 
therefore no need for simulating the half-way infinite lane segment leading up 
to that site. Allowing the integer I mS , if needed, to increase without bound, 
as we will do in the simulation, amounts to setting effectively L = oo. The 
reduced algorithm of the next subsection is based on these considerations. 

3.2 Reduced algorithm 

The reduced algorithm is an exact consequence of the full algorithm of section 
12.1.31 when the limit L — > oo is taken. This limit, almost paradoxically, 
simplifies the mathematics to the point that the reduced algorithm involves 
only the particle positions in the M x M intersection square and those in 
the row and column of entrance sites. There appear memory variables I mS 
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associated with these entrance sites. As the original injection sites have 
moved to minus infinity, particle injection now takes place de facto on the 
entrance sites. We state below only the points of difference between the 
reduced algorithm and the full one described in section 12.1.31 



3.2.1 Bulk step 

At the beginning of the sth time step, which covers the time interval s — 1 < 
t < s, the memory variables I m5 (s — 1) are known. The bulk step consists 
again of a sweep through all particles, ordered according to increasing phases, 
but now all located either on the M x M square or on one of the entrance 
sites. Particles move as in the full algorithm of subsection I2.1.1( the 'general 
case' and the 'first special case' of that section are the same here. 

Special case. The only special case to be discussed concerns the entrance 
sites. If during the sth time step the sweep encounters a particle j that is 
blocked on an entrance site (m, 5), then I mS (s) is calculated from I m5 (s — 1) 
according to (13. lap . 

If during the sth time step, say at time s — l + Tj, the sweep encounters a 
particle j that advances from the entrance site into the intersection square, 
then a random T™* is drawn as in the full algorithm and I mS (s) is deter- 
mined according to fl3.1bj) . The difference I mS (s — 1) — J m<5 (s) represents the 



reduction - compared to the free flow situation - of the number of time steps 
separating the departure of j and the arrival of j + 1 on the entrance site. 
Hence, recalling that s — l + Tj is the time at which j leaves the entrance 
site, we find that j + 1 must be injected onto that site at time tj + i given by 

t j+1 = s-l + T,+Tp + { + I mS (s)-I mS (s-l) 

= s-l + Tj + T^-mmil^is-l),^]), (3.2) 



where to pass to the second line we used ( 13.1bj) . Equation (13. 2p may be 



compared to (12. 2p . The time step of injection is s J+ i = L^j+iJ + 1 — s + L T j + 
T^_{\ - min (l mS (s - 1), L r J+iJ)- Particle j + 1, its phase r j+1 , its time step 
of injection Sj + i, and its lane index (m, 5) are placed on a list of particles 
waiting to be injected. 

Finally, an entrance sites that is empty is not involved in the sweep. For 
the corresponding lane (m, 5) the memory variable remains unchanged, that 
is, I mS (s) = I mS (s — 1), in agreement with (I3.1cp . 



3.2.2 Injection step 

Particle injection now takes place on the 2M entrance sites, to each of which 
the injection procedure is applied independently. If for a given entrance site 
(m, 5) the index s of the current time step is equal to the time step index 
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of the next particle to be injected on that site, say k, then the injection is 
carried out; the phase is inserted in the update sequence; and the injection 
is considered to have occurred at time s — 1 + Tfc. 

3.2.3 Initial state 

The initial state is generated as for the full model, section 12.1.31 but with L 
replaced by 1, that is, the injection sites coincide with the entrance sites and 
step (3) of section [2. 1.31 is empty. We set I rn5 (0) = for all (m, 5). 

4 Jamming transitions 

Whereas the preceding sections have dealt with the microscopic algorithm, 
the present section is of a theoretical nature. We derive certain relations that 
connect the averages directly obtained in the simulation to other physically 
meaningful averages. As observed in the remark that closes section 12.1.31 
both the initial state and the dynamics are determined by the 2M sequences 
of independent interval variables J™ 5 . An average is therefore a mean value 
calculated or measured with respect to the set {T™ 6 }. In a state with sta- 
tionary currents^) an average over a sufficiently long period of time must be 
expected to coincide with this J™ 5 average. Of primary interest will be the 
average currents Jjjf (a) for lane index m — 1, 2, . . . , M and as a function of 
the injection probability a. 

4.1 Current and reflection coefficient 

For each lane index m we must envisage two possibilities whose actual oc- 
currence is to be ascertained by the simulation. 

(a) The incoming flow is weak enough so as to pass entirely, which means 
that J,Jf (a) = J (a). The lane is then in a free flow phase. 

(b) The intersection square cannot handle the full incoming flow, which 
means that J^f (a) < J frcc (a). The lane is then in a jammed flow phase. 

We will write for either case 

= (i - R^V free , (4.i) 

where R^f (a) is the reflection coefficient of the mth lane and has the proper- 
ties (11. 3p . The product R^f J may be interpreted as the reflected current in 
the mth lane, to be discussed in greater detail in section E3~3l Since J free (a), 
given by (]2.4p . is imposed and the J^f (a) may be measured by counting the 

6 To which we will refer as the 'stationary state', even though the memory variables, 
and concomitantly the waiting line, may or may not be stationary. 
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outgoing particles in each lane, substitution of these two quantities in (14. ip 
yields the R* f . 

There is, however, a direct way of determining the R^f in the simulation. 
We will show below that in the long time limit the stochastic memory variable 
I m5 (s) satisfied 

I rn5 (t) 

l im _A! = R^ m = l,2,...,M, (4.2) 

t— YOO t 

so that R^f is also the average rate of growth of this variable. 



4.2 Equivalence of (HID) and (EES) 

We will take (14. 2p as the definition of the R^f and show that (14. ip follows. In 
the case of free flow the waiting line fluctuates only within a finite localization 
length, hence I mS (s) remains effectively bounded whence R^f = 0, and (14.11) 
is trivially true. It suffices therefore to consider the case of a jammed phase, 
for which I mS (t) is asymptotically linear in t. There are three steps to the 
proof. 

First step. In every time step s one of the three equations (13.11) is applied 
to obtain I mS (s) from I mS (s — 1). In the stationary state, let fs, f^, and 
/e be the fractional of all time steps in which equations (13. lap . (13.1bp . and 
(I3.1cp . respectively, are applied. When I mS (t) grows without bounds, the 
maximum m equation fl3~Tb|) is equal to I mS {s - 1) - [T™{\. It follows that 
for large times t we have the asymptotic proportionality 

I m5 (t) ~ {Oxf E + lxf B -KXf A )t, 

= (f B -Kf A )t, (4.3) 

in which n is the average of an d is easily calculated as 

f°° 1 - a 
k= / dT[T\P(T) = (4.4) 

Jo « 
with P{T) given by (123]) . Using (143]) . and (03]) we deduce that 

R„ A f = f B -—fA, (4-5) 
a 

which completes the first step. 

Second step. Using results of earlier work [25] we now determine the event 
fractions /a and fs in terms of a and /3^f. We consider a particle on the 

7 Almost surely, in the mathematical sense. 

8 The notation leaves the dependence of these fractions on m and M implicit. 
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entrance site. Following [25] we let a parameter 1 / /3^f stand for the average 
number of time steps that such a particle has to wait before it can enter the 
intersection squaraj. A particle that leaves the entrance site during the sth 
time step may or may not be replaced during the same time step. It was 
shown in reference [2S] that if, in the jammed phase, this site is unoccupied 
at some integer time s, then it will certainly be reoccupied at time s + 1. Let 
v be the average number of particles that cross the entrance site^ between 
two successive integer instants of time at which this site is unoccupied. This 
number v is determined exclusively by a and given by [25] 

1 = 1 + 1-1. (4.6) 
v a a 

It follows that for each time step at which the entrance site is unoccupied 
(event E above), there are on average v/(3^ time steps at which it is occupied 
(events A and B above). Out of the latter, there are v time steps at which 
the particle advances (event A), and a remaining (l//3^f — 1J v time steps at 
which it stays blocked (event B). Hence when a lane is in the jammed state, 
the event fractions /a, and /# are given by 

f - 1 f - V f - MPS - 1 > U7) 

Substituting ( 14. 71) in (14.5p we obtain 



DM = V Pm ( 1 _ I \ (AO\ 



Only for positive R^f is this result consistent with our initial supposition that 
lane m is in the jammed phase. Therefore, (14. 8 \j shows that lane m is jammed 
when a > (3^. Remembering that (3^ is a function, although unknown, of a 
and M we see that the jamming point a = is the solution of (3^ (a) = a. 



Third step. In the jammed phase the particle density on any site of the 
waiting line, and in particular on the entrance site, is given by f^ am = Ja + Jb 
and the current entering the intersection square is equal to /3* / / o' am . Since in 
the stationary state, the current entering the intersection square in lane m 
is equal to the outgoing current J^f in that lane, it follows with the aid of 

9 The inverse /3* f is the probability per time step that a particle on the entrance site is 
allowed to enter the intersection square. In reference [25] this parameter (called /? there) 
was an independent control parameter and entrance events of successive particles were 
uncorrelated. In the present case /3^f is a complicated function of a and M; moreover, 
correlations must be expected between entrance events. 
10 This sequence of particles is called a 'platoon'. 
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Figure 3: Lane containing a domain wall that separates a free flow from a jammed 
flow region. The particles in the latter constitute the waiting line. The particle on 
the entrance site (shaded) is the reference particle for the lane's memory variable. 
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(4.9) 



which is the jammed phase counterpart to (12. 4ft . Upon eliminating /3^f from 
(14. 8 p and (14. 9ft . solving for Jjjf in terms of R^f, and using (12. 4p . we obtain 
(14. ip . We have therefore established the interpretation of R^f, initially de- 
fined by (14. 2p . as the reflection coefficient of the mth lane. This completes 
the equivalence proof of this subsection. 

Equations (I4.ip . (14. 8p . and (14. 9 p show that out of the three quantities 
R^f , J^ f , and /3^f , each one determines the two others. 



4.3 Reflected current and domain wall motion 

We consider a single lane m and suppose that a > a^f so that jamming will 
occur. We now refer to figure |3j Going from left to right, the first particle 
in the waiting line is by convention the leftmost particle ever to have been 
blocked; in the figure this is particle j — 10. We will say that a 'domain wall' 
is located immediately to the left of this particle. The waiting line constitutes 
a finite spatial domain of jammed flow, separated by the domain wall from 
a half-infinite region in the free flow phase. The domain wall motion in a 
finite one-dimensional system has been discussed in references [25| [28]. 

Our reduced algorithm no longer contains the description of the domain 
walls. We will show below, however, that if at time s we know I mS (s) for some 
lane (m, 8), we can relate this quantity analytically to the average position 
that the domain wall in that lane would have at that time. For the initial 

11 This cannot be concluded from the figure alone, since the property of having never 
been blocked depends on the history of the configuration. 
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conditions of section 13.2.31 the domain wall is located at the entrance of the 
intersection square and, for a > a^f, starts at time t = propagating in 
the negative direction at some average speed that we will calo v^f . If we 
assimilate the incoming particle flow to a wave, then the domain wall is the 
moving front of the reflected wave. 

Consider now a time T large enough for fluctuating variables to be ap- 
proximated by their averages. Let L^f be the linear size of the jammed flow 
region, so that v^f = L^/T. To determine we reason as follows. The 
average number of particles N^f that up to time T has been prevented to 
cross the interaction square is 

N M = T(jfr ec_ J M ) _ (410) 

These particles are spread out along the waiting line and therefore the density 
^jam Q j? jj ne j g e q Ua ] pfree pj ug extra contribution N^f/L^f. This 

leads to a continuity equation in the form 

^ ai >, ft) = P fre » + t~m [</ fre » - O] , (4-11) 

m 

in which p> am (a, /3* f ) and p frec (a) are known from (14.91) and (12.41) . respectively. 
Using that v^f = L^/T we deduce that the speed of propagation v^f of the 
domain wall is given by 



v 



„free QM Jam 

M r ' rn r 

pjam pfree 

RM\ (l _ QM 



- - 1 . (4.12) 
a J \ v a J 

The second line results from some algebra in which (14.91) and (12.41) are used; 
the expression may be rewritten in several other ways. When eliminating /3^f 
from (I4.12p and (14.81) we obtain the speed of propagation v^f of the reflected 
wave as a function of R,^f. Explicitly, 

v M = ^> , (4.13) 

in which the coefficients are known functions of only the particle injection rate 
a. This equation shows that we may retrieve the average domain wall position 
from the reduced algorithm, in spite of the fact that at the microscopic level 
the waiting line has been eliminated from the description. Equation (14. 13ft 
is valid only for a > a^f . When a decreases to , we have that R^f —> 
and hence ~ [a/ (1 - a))R^ 0. 



■Quantities related to this wall will be indicated by sans serif symbols. 
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Figure 4: Reflection coefficient R\(a) for M = 1. The Monte Carlo algorithm 
of this work, which simulates an infinite system, agrees very well with the known 
analytic result [2D]. 

5 Monte Carlo results 

Using the reduced algorithm of section [3] we have carried out simulations 
of the intersecting streets for values of M up to M = 24. We determined 
the phase diagram of the stationary state of the intersection square as a 
function of the injection probability a and the linear size M. In our figures 
we will present the reflection coefficients R^f(a), which by (11. 2p are directly 
equivalent to the currents J^(a). Each simulation was started at t — 
from the initial random free flow configuration described in section 13.2.31 
Relaxation to a stationary state appeared to be very rapid. 

It is important to stress that, since we are in the limit L = oo, the 
simulation is free of finite size effects. The remaining statistical errors in the 
simulation results are entirely due to the finiteness of the simulation time. 

5.1 Phase diagram for M = 1 and M = 10 

For M = 1 the intersection square consists of a single site and there is a 
single critical jamming point a\ = | as predicted in earlier work [20]. The 
M = 1 simulation involves less than a dozen variables, in spite of the fact 
that we simulate an infinite system. Figure H] shows the simulation results for 
the reflection coefficient R^a). Each data point results from an average over 
1.1 x 10 6 time steps and over the x and y direction. One observes that R} is 
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Figure 5: Reflection coefficient R^(a) for M = 10 and m = 1,2,..., 10. The 
lines connecting the points are meant to guide the eye. The thin dotted vertical 
lines correspond to the critical points on the a axis; in increasing order these are 
ch\q, ctg , . . . , a\°. The thin dotted horizontal lines are the approximate levels at 
which the curves have their first and second discontinuity of slope. 

nonzero only for a > ~ and that it tends to | with infinite slope when a — > 1. 
A different representation of figure H] was obtained analyticalhl^l by Appert- 
Rolland et al. [20], whose results J free = a/(l + a) and J\ = v/(2v + l) imply 
that R{ = (u — u/a + 1)/ {2v + 1). This theoretical curve is also shown in the 
figure; its full agreement with the Monte Carlo data confirms the correctness 
of the 'reduced' algorithm of this paper. 

Figure [5] is for M — 10 and presents the reflection coefficients R™(a) for 
m — 1, 2, . . . , 10. Each data point results from an average over 1.1 x 10 7 time 
steps and the statistical error is less than the symbol size. There are ten 
critical values < al° < . . . < a\°. For a^ +l < a < a]^ the outer lanes 
1,2, ... ,m are in a free flow phase and the inner ones m + 1, . . . , M — 1, M 
have jammed flow. For each m the data shown are averages on the two 
lanes (m, x) and (m, y). We verified that no spontaneous symmetry breaking 

13 The M = 1 system actually studied in [20] was more general: it had four control 
parameters (two entrance and two exit rates) which could break the symmetry between 
the two crossing lanes; in the present work the two streets are symmetric. 
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Figure 6: Inverse of the principal critical point Cfe^j clS cl function of the linear size 
M of the intersection square (notice the logarithmic scale of M). The errors in 
the data points are smaller than the symbol sizes. The straight line is the best 
asymptotic linear fit. 

occurs between the x and the y direction. 
5.2 Phase diagram for general M 

We performed simulations for M = 1, 2, . . . , 10, 12, 15, 20, 24. Our findings 
for the other M values are qualitatively similar to those obtained in figure 
[5] for M = 10. Again, there is no spontaneous symmetry breaking between 
the two perpendicular directions. As a increases, the system first reaches 
its principal critical point at which the two innermost lanes (with m = 
M) get jammed. Then there is an initially fast succession of critical points 
< < a M-2 < • • • at each of which a further pair of lanes, one in 

each street, gets jammed. The spacing between the critical points becomes 

the outermost pair 



M 

i • 



gradually larger and at the last critical point, a = a 
of lanes (with m — 1) gets jammed. 

The principal critical point decreases with M but its exact law is 
difficult to ascertain. It is very well approximated by 



a M — 



A + B log M 



M 



> 



(5.1) 



as shown in figure El where A = 1.287 and B = 2.306. The high precision of 
these results is due to the fact that even if the street widths M are finite, the 
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simulations correspond to infinite street length, L = oo, and hence there are 
no finite length effects. The question of the large M limit of the critical point 
was asked also in the context of the BHL model [2T1 IM] . All these authors 
have declared to be unable to state whether or not this point goes to zero 
in the limit M — > oo. On the basis of our above results, albeit for relatively 
small M values, one might guess that (15 .ip is the correct asymptotic law and 
hence that for the model of this work the principal critical point does tend to 
zero with increasing M. We refrain at present, however, from drawing this 
conclusion. The reason is that preliminary simulations for larger M reveal 
complications in the form of metastabilities associated with the transition 
possibly turning first order. Since the infinite system limit is a fundamental 
question in statistical mechanics, we believe it is worthwhile, and necessary, 
to spend further efforts on studying it carefully. 

The intervals between two successive critical points scale in fairly good 
approximation as 

Ci - a m - °m + — , 5.2 

m z 

where Cm is a slightly negative and D M a positive constant. Equation (15. 2p 
appears to hold for m not too small and up to m = M. For M = 20 this is 
shown in figure [TJ the fit has C 2 o = —0.000 26 and _D 2 o = 0.605. 

For a > a^f the curve R^f (a) is composed of m segments that join with 
a discontinuity of slope at the higher critical values a^f-i; a m-2i ■ ■ ■ 5 F° r 
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M and m not too small, these segments, although actually curved, are very 
close to linear. The slope of the first segment of R^f(a), evaluated at ct^f, is 
approximately proportional to m 2 . Therefore, the larger M, the steeper the 
initial rise, and the finer should be the grid of points on the a axis in order to 
obtain a resolution of the successive jamming transitions in the inner lanes. 
This renders the determination of R™(a) by simulation increasingly harder 
as M grows. 

The first discontinuity of slope in each of the R^f (a) appears when R^jf (a) 
reaches a level r± that is, again in fairly good approximation, m independent. 
This had to be expected from the scaling observed above for the first segment 
of R^(ct), whose horizontal and vertical extension are ~ m~ 2 and ~ m 2 ), 
respectively. For M = 10 figure |5] shows that this level value is t\ ~ 0.47. 
Similarly, the second discontinuity appears at R^(cx) ~ r 2 ~ 0.69, the third 
one (not explicitly indicated in the figure) at r 3 « 0.83, and so on. The 
values of these levels appear to depend only little on M, as anticipated in 
the notation. 

As M grows, the values a^f with fixed m seem to tend to a limit value that 
may be called a^, but that we have not tried to determine with any precision 
in this study. Finally, for fixed a and M, the reflection coefficient approaches 
unity, and hence the transmitted current vanishes, roughly exponentially 
with m. It is still a challenge to find an analytic explanation for all these 
qualitative and quantitative observations, which were not expected a priori. 
Whereas an exact solution seems beyond reach, we think that an approximate 
theory may be possible. 

5.3 Snapshots 

Figure[8]is a snapshot of a 10 x 10 intersection square simulated with an injec- 
tion probability a = 0.169. It is the only one of our simulations that was car- 
ried out with the 'full' algorithm. The length of the incoming streets is equal 
to L — 15. For this value of a figure shows that the lanes m — 1,2, ... ,7 
are in the free flow phase whereas those with m = 8, 9, 10 are in the jammed 
flow phase. This is roughly what is visible in figure |HJ nevertheless, consid- 
erable fluctuations occur between independent snapshots, with lanes getting 
jammed and opening up intermittently. These fluctuations are enhanced by 
the fact that for finite L there are no sharply defined transition points, L 
being automatically an upper limit on the length of the waiting lines. 

Figure M is a snapshot of a 20 x 20 intersection square together with its 
row and column of entrance sites. The injection probability is a = 0.15. For 
M = 20 a phase diagram analogous to figure [5] shows that af± < 0.15 < a 2 [}. 
As a consequence, for a = 0.15 lanes 1 through 10 are in the free flow phase 
and 11 through 20 are in the jammed flow phase. 

Snapshots of this system taken at sufficiently long time intervals show 
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Figure 8: Snapshot of an intersection square of linear size M = 10 and for injection 
probability a = 0.169, obtained by the 'full' algorithm with L = 15. For this value 
of a, lanes 1 through 7 are in the free flow and 8, 9, 10 in the jammed flow phase. 
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Figure 9: Snapshot of an intersection square of linear size M = 20 for injection 
probability a = 0.15, obtained by the 'reduced' algorithm (L = oo). The row and 
column of entrance sites are also shown. For this value of a lanes 1 through 10 are 
in the free flow and 11 through 20 in the jammed flow phase. 
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considerable variation, but all have in common the existence of an approx- 
imately square (or sometimes rectangular) high density region in the lower 
left corner of the intersection square. The shape of this corner region fluctu- 
ates, since the jammed phase current still lets particles pass relatively easily 
along theb outermost jammed lanes (with indices m = 11, 12, ...) of this re- 
gion. However, the jammed phase current in the innermost lanes (those with 
indices m — 18, 19, 20) is very close to zero and the particles in those lanes 
are in a quasi-permanent frozen state. The net result is that there occurs a 
"freezing out" of a set of inner lanes, which reduces the traffic problem on 
the 20 x 20 intersection square to an effective one on a smaller square. In 
figure El this effective intersection square would correspond to a square of 
size 10 x 10 in the upper right corner. This reduction is certainly not exact 
but may well offer the starting point for a first theoretical approach. 

5.4 Further comments 

Our algorithm for crossing streets of infinite length was reduced to an algo- 
rithm on the finite-size intersection square with special memory boundary 
conditions. It is therefore natural that we make a comparison with existing 
models on finite lattices. A prominent one is the BHL model due to Biham 
et al. [21]. In this model two types of particles move unidirectionally on 
a torus of size M x M, one type horizontally and the other one vertically. 
In the original version of this system the horizontally and vertically moving 
particles are updated in parallel at the even and odd time steps, respectively. 
A single phase transition was observed between a free flow phase and a fully 
jammed (zero flow) phase. The importance of the aspect ratio of the lattice 
was stressed by D'Souza [29], who also showed that there exists a third, high 
density, jammed phase. In reference [21] and in the work that it has sparked 
[301 ED E21 E31 [29] , several variants of this BHL model were studied, many of 
them introducing additional stochastic elements. Of particular interest in our 
context is a recent brief report by Ding et al. [34] . These authors considered 
the BHL model with the standard open boundary conditions that consist in 
filling an empty site on the left or bottom boundary with a fixed probabil- 
ity at each time step. Under such boundary conditions there is no place for 
waiting lines and no sharply defined phase transition can be expected as long 
as M is finite. Obviously, the two types of boundary conditions correspond 
to distinct driving parameters and have distinct sets of applications: urban 
road traffic for the BHL model and intersecting pedestrian traffic flows for 
ours. 

The model of reference [31] also differs from ours by its use of random 
sequential update; but although the update type is important for the inter- 
pretation of the model, we consider that difference as secondary. 
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The most striking similarity between the results of reference [31] and 
ours is that when the driving parameter increases, the M x M domain is 
gradually filled with a dense phase that occupies a square or rectangular 
lower left corner, in the way shown in our figure M We find, as in reference 
[33] . that the complementary upper right corner, which is in a free flow state, 
has a size roughly independent of M. 

As for the differences, our M x M domain is part of two infinite streets 
and the memory boundary conditions keep track of the waiting lines in each 
of the 2M lanes. These boundary conditions create long-range correlations 
in time, that in turn determine sharply defined phase transition points on 
the a axis for all M = 1, 2, 3, . . .. As a consequence, in the model of this 
paper the progressive growth of the dense corner occurs through a sequence 
of M phase transitions. Our simulations concern relatively modest values of 
M and are aimed at locating the transition points with high precision. In 
addition we are able to assign to each of the 2M lanes a reflection coefficient 
which serves as a lane order parameter, and a waiting line of which at each 
instant of time the average length is known. 

6 Conclusion 

We have introduced and studied a lattice model of pedestrian traffic on two 
crossing one-way streets. Each street is represented by a set of M paral- 
lel TASEP lanes and the only model parameter besides M is the injection 
probability a of a pedestrian ('particle') at minus infinity. The dynamics 
is based on frozen shuffle update [23J [23 particles enter stochastically 
but once in the system move deterministically. From an algorithmic point 
of view we have found that in this model the frozen shuffle update leads to 
accurate results even with a modest simulation effort. We consider this as 
an encouragement for simulating other models, possibly very different ones, 
with the same update. 

The intersecting streets that we considered are infinite in both directions, 
but one achievement of this work has been to show that the dynamics may be 
reduced to a problem of interacting variables on the sites of the finite M x 
M intersection square. Appropriate boundary conditions were formulated 
in terms of 'memory variables' and we established the theoretical relation 
between these new quantities and the outgoing current. 

Our Monte Carlo work shows that as a increases, the system undergoes a 
sequence of M phase transitions starting with the one at the principal critical 
point a^. Since we are able to perform the simulation directly on an infinite 
system, there are no finite size effects; the uncertainty in the critical points is 
due only to the finiteness of the simulation time. At each transition a new pair 
of lanes, one in each street, passes from a free flow to a jammed flow phase. 
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We find that a reflection coefficient for the current is an appropriate order 
parameter. The accuracy of the algorithm has allowed us to establish certain 
surprising features, such as the discontinuities of slope that occur when the 
reflection coefficient crosses a set of narrowly defined level values. All these 
results should help trigger interest in building theoretical approaches to this 
and similar systems. 

The present study stays far from exploring, let alone answering, all ques- 
tions that this specific model poses. Quantities of interest not studied here 
include the particle density in the intersection square, its space dependence, 
and its fluctuations and correlations; the fluctuation of the waiting line 
lengths at and near the critical points; and its relation to the fluctuation 
of the memory variables. An unanswered question also concerns the behav- 
ior of the principal critical point as M — > oo. 

We believe that this model, because of its simplicity, sets a standard sce- 
nario with respect to which others may be discussed. Clearly, ideas come 
readily to mind about how to modify or extend this model, for example by 
considering streets of different widths Mi and M2, unequal injection proba- 
bilities a x and a y , and so on; or by opening the possibility for pedestrians to 
move laterally or diagonally forward. Sideways steps between the lanes are 
certainly expected to blur the distinction between the jamming transitions in 
the individual lanes that here occur separately. It is hard to guess, however, 
what sort of a transition will then result. Several effects of the present model 
may turn out to be robust, such as the predominance of the flow through the 
outer lanes over those through the inner ones. In future work [35] we will 
address a few of the many new issues raised here. 
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